Muscle coordination retraining inspired by musculoskeletal simulations reduces knee contact force

Humans typically coordinate their muscles to meet movement objectives like minimizing energy expenditure. In the presence of pathology, new objectives gain importance, like reducing loading in an osteoarthritic joint, but people often do not change their muscle coordination patterns to meet these new objectives. Here we use musculoskeletal simulations to identify simple changes in coordination that can be taught using electromyographic biofeedback, achieving the therapeutic goal of reducing joint loading. Our simulations predicted that changing the relative activation of two redundant ankle plantarflexor muscles—the gastrocnemius and soleus—could reduce knee contact force during walking, but it was unclear whether humans could re-coordinate redundant muscles during a complex task like walking. Our experiments showed that after a single session of walking with biofeedback of summary measures of plantarflexor muscle activation, healthy individuals reduced the ratio of gastrocnemius-to-soleus muscle activation by 25 ± 15% (p = 0.004, paired t test, n = 10). Participants who walked with this “gastrocnemius avoidance” gait pattern reduced late-stance knee contact force by 12 ± 12% (p = 0.029, paired t test, n = 8). Simulation-informed coordination retraining could be a promising treatment for knee osteoarthritis and a powerful tool for optimizing coordination for a variety of rehabilitation and performance applications.

between the medial and lateral compartments of the knee. However, these interventions do not target, and often increase, the large muscle contribution to KCF 19 . Additionally, despite the relationship between knee loading and pain, individuals with osteoarthritis do not naturally select coordination strategies that minimize joint loading [20][21][22] . Even when provided with real-time feedback of KCF, individuals are unable to alter their muscle coordination to reduce KCF 23 . This suggests that the complex dynamics relating muscle coordination to KCF, and not exclusively a lack of feedback, could limit the ability of individuals with osteoarthritis to adopt a loadreducing coordination strategy.
Musculoskeletal simulations can identify coordination strategies that reduce KCF. The second, or late-stance, peak of knee contact force (KCF P2 ) is sensitive to changes in muscle coordination and can, in theory, be altered by several times bodyweight without changing kinematics 21,24,25 . A coordination strategy that minimizes KCF P2 can be achieved by altering the activation of every lower-extremity muscle 21,25 , but this muscle coordination pattern is too complicated to learn. Solutions that involve simpler changes in muscle coordination, like reducing the activation of the gastrocnemius 12,21,[24][25][26] , are likely easier to learn and can still substantially reduce KCF P2 . While simulations suggest that a "gastrocnemius avoidance" coordination strategy could effectively reduce KCF P2 , it remains unclear which redundant muscles would need to compensate for the reduction in gastrocnemius activation. Furthermore, it is unknown whether humans can learn to change the relative activation of redundant muscles during a complex task like walking.
Electromyography (EMG) biofeedback is an effective tool for exploring the limits of volitional motor control. During walking, biofeedback of the tibialis anterior or gastrocnemius can help restore normal ankle dorsiflexion or plantarflexion moments for individuals with stroke or cerebral palsy 27,28 . Humans can also gain selective control over several motor units within a single biofeedback session 29,30 , but most individuals lose this selective control during dynamic tasks 31 . Additionally, EMG biofeedback can aid in selecting kinematic strategies that change the relative activation of redundant muscles during physical-therapy exercises 32,33 . There has been limited work, however, exploring whether humans can change the relative activation of redundant muscles without changing kinematics during dynamic tasks like walking.
The purpose of this study was to design and implement a muscle coordination retraining intervention that teaches individuals to reduce knee contact force. We first used musculoskeletal simulations to identify a subset of muscles to target with biofeedback that would reduce knee contact force without changing joint kinetics (Fig. 2a). Based on the simulation results, we designed an EMG biofeedback intervention to test whether healthy individuals could learn to change the coordination of redundant muscles during walking (Fig. 2b). We then used EMG-informed simulations to evaluate whether participants reduced their knee contact force when walking with the new coordination pattern (Fig. 2c). More generally, we describe a framework for using simulations to design simple biofeedback interventions that teach individuals new coordination patterns that have a therapeutic or performance-enhancing effect. Contributions to knee contact force during walking. Muscle force contributions to knee contact force (red) exceed the intersegmental force contribution (blue), which is the reaction force from inverse dynamics. Muscle force contributions are dominated by the hamstrings during the first 10% of stance, the quadriceps from 10-40% stance, and the gastrocnemius from 40 to 90% stance. Forces are shown in terms of bodyweight (BW) for the 10 healthy subjects in this study, and muscle activations were estimated by minimizing the sum of squared activations using static optimization. www.nature.com/scientificreports/

Results
Simulation-based intervention design. To design the muscle coordination retraining intervention, we sought to understand the compensatory muscle activations that are necessary to generate normal walking kinetics with minimal gastrocnemius activation (Fig. 2a). We simulated healthy walking kinetics with two different static optimization objective functions: a natural coordination objective that is a surrogate for metabolic energy 34,35 (Eq. 1) and a gastrocnemius avoidance objective that additionally penalizes gastrocnemius activation 21 (Eq. 2).
These simulations showed that walking without activating the gastrocnemius requires increased force generation from the soleus, hamstrings, and hip flexors (Fig. 3). Increased soleus force supplements the ankle plantarflexion moment generated by the gastrocnemius, but unlike the gastrocnemius, the soleus does not cross the knee, allowing it to plantarflex the ankle without compressing the knee. Increased hamstrings force supplements the knee flexion moment generated by the gastrocnemius, but the hamstrings have, on average, a knee flexion moment arm that is 1.7 times greater than the gastrocnemius 36 , allowing these muscles to generate the same moment with less force. The antagonistic hip extension moment generated by the hamstrings required an increase in iliopsoas force. The large changes in the force required by the ankle plantarflexors coupled with their essential role in generating normal walking kinematics 37,38 led us to focus the biofeedback intervention on the relative activation of the gastrocnemius and soleus muscles. With the goal of teaching participants to reduce gastrocnemius force without reducing their plantarflexion moment, we designed biofeedback that instructed participants to reduce gastrocnemius activation and increase soleus activation. Figure 2. Simulation-inspired biofeedback design process. (a) To design the biofeedback for teaching individuals to walk without activating their gastrocnemius (gastroc), we simulated walking using both a natural and a gastrocnemius avoidance objective function (Eqs. 1 and 2) in the static optimization muscle redundancy solver. (b) Based on the simulation results, we provided individuals with visual biofeedback of their muscle activation, measured with electromyography (EMG), that instructed them to change the coordination of their ankle plantarflexor muscles. Participants performed five walking trials: a baseline trial (natural walking), three trials with visual biofeedback, and a retention trial without feedback. During all three feedback (FB) trials, visual biofeedback instructed participants to reduce their gastrocnemius-to-soleus activation ratio (bar magnitude). During the final two feedback sessions, additional feedback was provided, instructing participants to also reduce their average gastrocnemius (gastroc) activity (bar color). (c) EMG-informed static optimization simulations were used to assess the effect of a gastrocnemius avoidance coordination pattern on knee contact force. The simulated gastrocnemius-to-soleus activation ratio was constrained to match the ratio measured with EMG. www.nature.com/scientificreports/ Muscle coordination retraining experiment. We conducted an experiment to explore whether humans could learn to change the relative activation of their gastrocnemius and soleus muscles during walking when given visual biofeedback. During a single session, 10 healthy adults performed five six-minute walking trials on an instrumented treadmill: a baseline trial, three feedback trials, and a retention trial (Fig. 2b). During all feedback trials, a real-time bar plot instructed participants to reduce their gastrocnemius-to-soleus activation ratio (Eq. 3): where EMG gastroc (t) and EMG soleus (t) are the average EMG linear envelopes over the stance phase for the medial gastrocnemius and soleus, respectively. During the final two feedback trials, participants were also given feedback to reduce their average gastrocnemius activation to prevent them from reducing the activation ratio by only increasing soleus activation. During the retention trial, participants were instructed to retain their learned coordination pattern but were not given feedback. After receiving visual biofeedback of their muscle activation ratio during the first trial, all participants reduced their activation ratio, by an average of 22 ± 12% (p < 0.001, t test, n = 10), but they did not significantly reduce gastrocnemius activation (p = 0.681, t test, n = 10, Fig. 4a). By adding gastrocnemius activation feedback during the second and third feedback trials, we sought to teach participants to reduce their activation ratio by reducing gastrocnemius activation. Following the third feedback trial, participants reduced gastrocnemius activation by 17 ± 19% (p = 0.033, t test, n = 10) compared to baseline. Finally, we investigated whether individuals could retain their new coordination pattern after six minutes of walking without feedback. At the end of the retention trial, participants retained a 25 ± 15% (p = 0.004, t test, n = 10) reduction in activation ratio and a 17 ± 21% (p = 0.033, t test, n = 10) reduction in gastrocnemius activation. Despite this change in activation ratio, ankle kinematics and kinetics remained similar to baseline (Fig. 4b-e). The average ankle moment only changed by 3 ± 14% during the retention trial compared to baseline, which trended towards statistical equivalence within one baseline standard deviation (p = 0.063, two-one-sided t tests for equivalence, n = 10).
EMG-informed simulation of knee contact force. We evaluated whether the gastrocnemius avoidance coordination pattern that participants adopted achieved the desired therapeutic objective: a reduction in simulation-estimated KCF P2 . Using static optimization, we simulated five gait cycles from the final minute of both the baseline and retention trials. To capture changes in ankle plantarflexor muscle activity between conditions, we constrained the gastrocnemius-to-soleus activation ratio at each time step in the simulation to match the corresponding ratio measured with EMG (Fig. 2c). The simulated activations qualitatively matched EMG patterns for the major muscle groups crossing the knee and ankle (Fig. 5).  Figure 3. Simulation of a gastrocnemius avoidance coordination pattern. Identical joint kinetics were simulated (n = 1) with a natural (Eq. 1) and a gastrocnemius (gastroc) avoidance (Eq. 2) static optimization objective function. During late stance, the muscles generate ankle plantarflexion, knee flexion, and hip flexion moments. The gastrocnemius avoidance coordination pattern requires increased soleus and hamstrings force to compensate for the ankle and knee moments normally generated by the gastrocnemius. The larger knee flexion moment arm (r) of the hamstrings compared to the gastrocnemius allows the hamstrings to generate the same moment with less force than the gastrocnemius, reducing knee contact force. A weighted average of the moment arms (hamstrings: biceps femoris long head and short head, semitendinosus, semimembranosus; gastrocnemius: medial and lateral heads) was computed at 75% of the stance phase, weighted by each muscle's optimal force in the musculoskeletal model 36 .  . Participants reduced their gastrocnemius-to-soleus activation ratio and average gastrocnemius (gastroc) activation during training. They retained these reductions following the retention trial (*p < 0.05, paired t test, p-values reported after controlling for the false detection rate). (b-e) The mean (line) and standard deviation (shading) of ankle plantarflexor muscle activity and ankle mechanics for the baseline (base.) and retention (ret.) trials. Despite the 25 ± 15% change in activation ratio from the baseline to retention trial, the stance-phase-averaged ankle moment only changed by 3 ± 14%, which trended towards being equivalent to baseline within one baseline standard deviation (p = 0.063, two-one-sided t tests for equivalence, n = 10).

Figure 5.
Simulated and measured muscle activation. Electromyography (EMG) linear envelopes averaged across all participants (n = 10) with a 40 ms electromechanical delay are compared to simulated muscle activations for the baseline and retention trials. Despite magnitude differences between EMG and simulated activations of the gastrocnemius and soleus, the relative changes between trials are consistent, indicating that our EMG-informed static optimization technique captured the changes in ankle plantarflexor muscle activity measured with EMG. The activations of the muscles in the top row were not informed by EMG in the simulation. The shape, timing, and between-trial changes in simulated activation matched EMG for these muscles, with the exception of the rectus femoris. www.nature.com/scientificreports/ To determine the mechanical efficacy of the intervention, we analyzed knee contact forces for the eight individuals who retained a reduction in late-stance gastrocnemius activation during the retention trial. On average, these individuals reduced their KCF P2 by 0.38 ± 0.39 BW (12 ± 12%, p = 0.029, t test, n = 8) during the retention trial compared to baseline (Fig. 6a). Six of the eight individuals who reduced late-stance gastrocnemius activation reduced their second peak knee contact force, but five increased their first peak (Fig. 6b). The 0.25 ± 0.31 BW reduction in the gastrocnemius contribution to KCF was the main contributor to reduced KCF P2 (Fig. 6c). Contrary to the prediction by the intervention-design simulation (Fig. 3), neither the hamstrings EMG linear envelope (Fig. 5) nor the hamstrings contribution to KCF (Fig. 6c) increased during late stance. The only significant difference in lower-extremity sagittal joint kinematics or kinetics was a 0.72 ± 0.63%BW*height (31 ± 24%, p = 0.015, t test, n = 8) reduction in the late-stance knee flexion moment (Fig. 7). This finding aligns with our observed reduction in gastrocnemius activation but unchanged hamstrings activation.

Discussion
We developed a muscle coordination retraining intervention that reduces knee loading. We first used musculoskeletal simulations to identify simple changes in coordination that reduce knee contact force, which informed the design of a biofeedback intervention. With simple biofeedback of their muscle activity, healthy individuals were able to quickly change the relative activation of their gastrocnemius and soleus muscles while still generating a similar ankle plantarflexion moment. This change in coordination resulted in a 12% reduction in simulationestimated knee contact force, on average, that may have therapeutic benefit for individuals with knee osteoarthritis. These results suggest that simulation-based intervention design is a promising tool for identifying new, learnable coordination strategies that achieve a therapeutic objective.
Our results build upon prior EMG biofeedback work by demonstrating that humans can quickly alter the coordination of redundant muscles during a complex task with simple feedback. Acute coordination changes have been learned within a single visit for static activities 29 , but clinical coordination retraining programs can last many weeks 32 . All 10 of our participants reduced their gastrocnemius-to-soleus activation ratio during walking after only five minutes of feedback. This suggests that multiple coordination strategies could be taught during a single visit for applications, like human-in-the-loop optimization, in which users need to rapidly adapt to new conditions 39 . Prior work has also demonstrated that precise control of individual motor units or muscles is more easily attained during isometric contractions 29,31 or controlled physical therapy activities 32,33 , but this study suggests that the relative activation of redundant muscles can be retrained during dynamic activities like walking. Further studies are necessary to determine how complex of coordination changes can be taught during functional activities.
Pizzolato et al. showed that during a single visit, healthy individuals were not able to identify coordination patterns that reduced medial knee contact force with real-time feedback of this force alone 40 . This could be due to the complexity of the dynamics that relate muscle coordination to KCF. When given feedback of the activation of two muscles, most of our participants reduced KCF, which highlights the importance of selecting biofeedback targets that are both simple and effective.
During a single visit, six of 10 healthy individuals adopted the gastrocnemius avoidance coordination pattern in a way that reduced their late-stance knee contact force. These promising results can inform future longer-term studies of this intervention in individuals with knee osteoarthritis. First, age and pathology-related changes in Figure 6. The effect of a gastrocnemius avoidance coordination pattern on knee contact force. (a) The mean (line) and standard deviation (shading) of knee contact force for the participants who reduced late-stance gastrocnemius activation during the retention trial (n = 8). These participants reduced their second peak of knee contact force compared to baseline (*p = 0.029, paired t test). (b) Changes in the first and second peak contact force between baseline (base.) and retention (ret.) are shown for all participants (n = 10), with the two participants who did not retain a reduction in late stance gastrocnemius activation represented with dashed lines. Six of the eight individuals who reduced late-stance gastrocnemius activation reduced their second peak knee contact force, but five increased their first peak. (c) For the eight subjects who reduced late-stance gastrocnemius activation, the change in knee contact force is decomposed into the intersegmental and muscle force components. Reductions in the second peak of knee contact force are primarily driven by reductions in gastrocnemius force. www.nature.com/scientificreports/ strength, stability, and motor control may change how easily individuals with osteoarthritis can learn this intervention and how effectively it reduces knee contact force. Multiple visits may be required to evaluate if individuals can alter their coordination pattern, and many retraining bouts will likely be necessary for long-term retention (e.g., 6-8 sessions have been effective for retaining a kinematic gait modification 18,41 ). The variance in the intervention's effect on knee contact force in healthy individuals suggests that future clinical efficacy studies should ensure that individuals reduce their contact force prior to adopting the intervention for a long duration 42,43 . Additionally, understanding the effect of this intervention on measures of joint loading beyond KCF P2 is important prior to studying this intervention in individuals with osteoarthritis. On average, the intervention did not increase peak hip or ankle contact forces in the leg that received feedback (see "Effect of gastrocnemius avoidance gait on hip and ankle contact forces" in Supplementary Information). In this leg, the intervention reduced the second peak of medial knee contact force on average but increased the first peak of lateral knee contact force (see "Effect of gastrocnemius avoidance gait on medio-lateral knee loading metrics" in Supplementary Information). Five of eight participants increased the first peak of KCF in the leg that received feedback (Fig. 6b), and on average, participants increased their early-stance knee flexion angle and knee extension moment peaks in the leg that did not receive feedback, potentially increasing knee contact force (see "Effect of gastrocnemius avoidance gait on the contralateral limb" in Supplementary Information). An increase in the early-stance knee extension moment or KCF peak could put individuals with osteoarthritis at risk of accelerated cartilage degeneration 9,44 . Providing bilateral electromyography biofeedback or biofeedback discouraging increases in the knee extension moment could attenuate these increases in early-stance knee loading. Future studies that investigate this coordination retraining intervention in individuals with osteoarthritis should test these strategies for mitigating elevated early-stance knee loading, and the intervention should likely not be prescribed to individuals who cannot adopt it without a substantial increase in early-stance knee loading in either leg. Simulation-inspired coordination retraining is an effective and flexible way to alter joint loading. The individuals who adopted the gastrocnemius avoidance coordination pattern reduced KCF P2 by an average of 0.38 BW, which is similar to the joint-offloading effect of losing 15-38% of bodyweight [45][46][47] . Similar reductions in medial and total KCF have been achieved with kinematic gait modifications 48,49 or assistive devices 15,48,50 , but coordination retraining does not require conspicuous kinematic changes or lifestyle changes that often hinder patient adoption 18,51 . More generally, coordination retraining reduces the muscle contribution to joint loading by leveraging functional differences commonly found between redundant muscles. This approach could target the total, medial, or lateral KCF, making it more flexible than current interventions that primarily change the medio-lateral distribution of KCF by reducing the knee adduction moment 14,15,[17][18][19] . Coordination retraining could even be coupled with kinematic modifications (e.g., toe-in gait) to reduce medial contact force during both early and late stance. Furthermore, coordination retraining could reduce loading in other joints, like the hip 25 , for which there are a limited number of offloading interventions.
In addition to designing interventions to reduce joint contact force, musculoskeletal simulations can identify coordination retraining interventions that optimize for other important clinical and sports-performance metrics. Simulations can identify optimal muscle coordination strategies for walking or running with a robotic assistive device 52,53 . Providing simulation-guided EMG biofeedback may augment the natural motor learning process, allowing a device user to more quickly converge on an energetically favorable coordination strategy. Additionally, Figure 7. Changes in knee and hip mechanics. The mean (line) and standard deviation (shading) of knee and hip kinematics and kinetics for the participants who retained a reduction in their late-stance gastrocnemius activation during the retention trial (n = 8). During the retention trial, participants walked with a smaller latestance knee flexion moment (*p = 0.015, paired t test). www.nature.com/scientificreports/ muscle coordination patterns identified by simulations could improve sports performance. For example, coordination patterns that maximize the distance of a long jump 54 could be taught to athletes with simplified EMG biofeedback. Finally, this framework has potential applications in injury prevention and rehabilitation. For example, simulations may be able to identify new coordination patterns that reduce strain in commonly injured ligaments, like the anterior cruciate ligament or ulnar collateral ligament (e.g., by changing the activation of ligament agonists and antagonists [55][56][57] ). Even after an injury, new muscle coordination patterns may help normalize joint mechanics 24 , potentially making coordination retraining an important part of post-injury rehabilitation. It is important to acknowledge the limitations of this study. We demonstrated that eight of 10 healthy individuals could learn and retain a new coordination pattern in the context of a single visit; further work is necessary to evaluate learning and retention over long durations in free-living conditions. Although there are challenges with day-to-day EMG normalization, EMG-based biofeedback is amenable to a wearable solution for in-home training. Second, we show the feasibility of retraining muscle coordination to reduce knee contact force in young, healthy participants, but it is unclear how effective this intervention would be for individuals with knee osteoarthritis. Future studies should investigate whether individuals with osteoarthritis can learn a gastrocnemius avoidance gait pattern and whether it reduces both knee contact force and pain. Third, our participants, on average, did not maintain their natural walking kinetics when adopting the gastrocnemius avoidance coordination pattern. Our feedback successfully taught individuals to change ankle plantarflexor muscle activity while maintaining a natural ankle plantarflexion moment; however, we did not give feedback to increase hamstrings activation, which may have resulted in the observed reduction in the late-stance knee flexion moment. If future studies aim to retain identical joint kinematics and kinetics, additional feedback should be given to instruct participants to achieve all of the necessary compensatory changes in muscle activation. Fourth, we estimated changes in KCF using a musculoskeletal model instead of directly measuring them with an instrumented knee implant. However, a similar static optimization approach has been shown to accurately detect gait modification-induced changes in KCF measured with an instrumented knee implant 58 . The changes in simulated muscle activation between natural and gastrocnemius avoidance gaits also matched changes in EMG (Fig. 5), providing confidence in the accuracy of our simulated changes in KCF. Finally, our intervention reduced the second peak of total and medial contact force on average, but many studies that relate knee loading to osteoarthritis progression have studied the first peak or impulse of the knee adduction moment 44,59 or the first peak of medial contact force 9 . Future studies should investigate which measures of knee loading are most related to osteoarthritis initiation and progression.
In summary, eight of 10 healthy individuals in our study learned and retained changes to their muscle coordination pattern during walking when given simple biofeedback. By altering the activation of their redundant plantarflexor muscles, most of these individuals reduced their late-stance knee contact force, demonstrating that coordination retraining may be a valuable new conservative treatment to offload osteoarthritic joints. More generally, simulation-inspired muscle coordination retraining may be an effective way to teach people novel, therapeutic coordination strategies. This approach capitalizes on the adaptability of the human neuromusculoskeletal system to optimize muscle coordination, with compelling injury prevention, rehabilitation, and sports performance applications.

Methods
Simulation-based intervention design. Previous work had suggested that knee contact force could be reduced with a gastrocnemius avoidance coordination pattern 21 , but it was unclear what compensatory changes in muscle activation were necessary. To address this question, we simulated one gait cycle of normal walking kinematics (healthy male, speed 1.25 m/s) with two feasible muscle coordination patterns: a natural pattern that minimized a surrogate for energy expenditure and a gastrocnemius avoidance pattern that minimally activated the gastrocnemius (Fig. 2a). We used OpenSim 4.0 60,61 and a custom static optimization implementation in MATLAB R2017b (Mathworks, Inc., Natick MA, USA) for musculoskeletal modeling and simulation.
We modeled the lower extremities and torso using the musculoskeletal model described by Rajagopal et al. 36 with 21 degrees of freedom, 40 musculotendon actuators for the lower extremity that received biofeedback, and nine ideal torque actuators for the contralateral limb and torso. The model included six degrees of freedom between the pelvis and the ground, three rotational degrees of freedom between the pelvis and torso, three rotational degrees of freedom at the hip, one rotational degree of freedom at the knee that parametrized the remaining rotational and translational degrees of freedom of the tibiofemoral 62 and patellofemoral 63 joints, and one rotational degree of freedom at each of the ankle and subtalar joints. We modified the model by calibrating each muscle's passive muscle force-length curves so that joint moments generated by passive muscle forces more closely matched experimental data 64 (see "Passive muscle force calibration" in Supplementary Information). We also altered the muscle paths of the hip abductor musculature to more closely match moment arms estimated from experiments 65,66 , finite element models 67 , and MRI 68 (see "Hip abductor muscle path adjustments" in Supplementary Information). The modified generic model was scaled using the OpenSim Scale Tool to match the anthropometric measurements taken from a standing static trial, and virtual markers on the model were moved to match experimental marker locations during this trial. Model scaling adjusts muscle optimal fiber lengths, tendon slack lengths, and muscle moment arms based on subject-specific, proportional segment scaling factors (the ratio of tendon slack length to optimal fiber length is preserved for each muscle). No other musculotendon parameters were adjusted. Model kinematics were estimated using the Inverse Kinematics tool in OpenSim, which minimizes the error between the positions of experimental markers and virtual model markers. Joint moments were computed using the Inverse Dynamics tool in OpenSim with low-pass-filtered kinematics (6 Hz, 6th order, zero-phase shift Butterworth in OpenSim) and ground reaction forces (6 Hz, 4th order, zero-phase shift Butterworth) as inputs. www.nature.com/scientificreports/ We developed a custom static optimization algorithm in MATLAB, using the OpenSim API, to solve the muscle redundancy problem. Two objective functions were used. The first minimized the sum of squared muscle activations (Eq. 1), which is a commonly used surrogate for metabolic cost 34,35 . The second was a gastrocnemius avoidance objective function with a heavily weighted penalty on gastrocnemius activation 21 (Eq. 2). Musclegenerated joint moments were constrained to match joint moments from inverse dynamics (Eq. 4) at each timestep ( t k ).
where M ID,j are the inverse dynamics joint moments for each of j degrees of freedom (DOF), r MTU i,j are the moment arms 69 of the ith musculotendon actuator (MTU) about the jth DOF, and F MTU i are the MTU forces. MTU forces were computed using the Hill-type model 70 described by Millard et al. 71 . We approximated tendon compliance by solving the musculotendon static equilibrium equation (Eq. 5) for muscle fiber length using the MTU length from the current step and the activation from the previous step: where F m o is the maximum isometric muscle force, f T l M , l MTU is the tendon force-length multiplier as a function of the muscle fiber and MTU lengths, a is muscle activation between 0 and 1, f l (l m ) is the active muscle fiber force-length multiplier, f PE (l m ) is the passive muscle fiber force-length multiplier, and α(l m ) is the muscle pennation angle. After solving Eq. (5) for l m (t k ) using Newton's method 71 , we fixed l m (t k ) and solved for F MTU (t k ) as a function of the design variable, a(t k ) , for each muscle (Eq. 6).
We then used fmincon in MATLAB to solve the static optimization problem with the combination of Eqs. (4) and (6) as a constraint and either Eq. (1) or Eq. (2) as the objective function. The static equilibrium computation in Eq. (5) does not account for the force-velocity property of muscle. However, in a pilot sensitivity analysis, we found that static optimization solutions that used a compliant tendon but excluded the force-velocity property of muscle yielded gastrocnemius and soleus activations that more closely matched EMG than solutions that assumed a rigid tendon but incorporated the force-velocity property (see "Tendon compliance sensitivity analysis" in Supplementary Information). This aligns with findings from dynamic simulations showing that Achilles tendon compliance enables the ankle plantarflexors to generate force with small shortening velocities during stance phase 72 .
Using the muscle forces from static optimization, we used the Joint Reaction Analysis tool in OpenSim to compute KCF as the reaction force along the longitudinal axis of the tibia. We also compared the differences in muscle forces in the major muscle groups that generate the late-stance ankle plantarflexion, knee flexion, and hip flexion moments: the soleus, the gastrocnemius (medial and lateral heads), the hamstrings (semitendinosus, semimembranosus, biceps femoris long and short head), the iliopsoas (iliacus, psoas), and the rectus femoris. Finally, after evaluating the changes in muscle force between coordination patterns, we chose to focus biofeedback on the two muscle groups with the greatest changes in force: the gastrocnemius and soleus. By reducing gastrocnemius activation and increasing soleus activation, an individual could, in theory, reduce the gastrocnemius contribution to knee contact force without changing the net ankle plantarflexion moment.
Muscle coordination retraining experiment. Eleven individuals enrolled in the study, and 10 (sex: 4 female, age: 26 ± 4 years, BMI: 22.8 ± 2.1) completed it. All individuals provided informed consent prior to beginning the experiment. The Stanford Institutional Review Board (IRB00000351) approved of the consent process and experimental protocol prior to participant enrollment. We performed the experiment in accordance with this approved protocol and all relevant guidelines and regulations. We included individuals who did not have a history of lower-limb injury in the past year and who did not have a history of knee ligament or meniscus injury. Before individuals began the walking portion of the study, we evaluated whether our surface EMG electrodes could measure different signals between the medial gastrocnemius and soleus. We asked participants to perform standing and seated plantarflexion exercises, since these activities induce large changes in the relative activation of the gastrocnemius and soleus muscles due to varying amounts of knee flexion 73 . Because our target reduction in gastrocnemius-to-soleus activation ratio during walking was 15%, we only included participants who reduced their activation ratio by at least 15% during the seated plantarflexion activity compared to the standing plantarflexion activity. One participant was excluded for this reason.
Participants completed a single visit to a motion capture laboratory with a force-instrumented treadmill (Bertec Corporation, Columbus OH, USA), an 11-camera motion capture system (Motion Analysis Corporation, Santa Rosa, CA, USA), and a wireless surface EMG system (Delsys Corp., Natick, MA, USA). We used custom MATLAB scripts for real-time computation and biofeedback. The scripts used force, marker, and EMG signals streamed from the motion capture software (Cortex v.7.0, Motion Analysis Corporation, Santa Rosa, CA, USA). We placed markers bilaterally on the 2nd and 5th metatarsal heads, calcanei, medial and lateral malleoli, medial and lateral femoral epicondyles, anterior and posterior superior iliac spines, acromion processes, sternoclavicular joints, and on the C7 vertebrae. Markers on the medial femoral epicondyles and malleoli were removed prior to walking trials, and 16 additional markers were used to aid in limb tracking. After randomly selecting a leg for analysis and biofeedback, we placed nine EMG electrodes unilaterally on the soleus (medial aspect), medial www.nature.com/scientificreports/ gastrocnemius, lateral gastrocnemius, tibialis anterior, biceps femoris, semitendinosus, rectus femoris, vastus medialis, and vastus lateralis. Prior to walking, participants performed maximum voluntary contraction activities for EMG normalization, static and dynamic calibration trials to aid in scaling of the musculoskeletal model, and calf-raises to evaluate EMG sensor placement. After walking on the treadmill to warm up, participants performed two resisted isometric and isokinetic maximum voluntary contraction activities: prone knee flexion and supine ankle dorsiflexion. Participants then performed five maximum effort hip flexion kicks and five maximum height jumps 74 . We calculated EMG linear envelopes by bandpass filtering (30-500 Hz, 4th order, zero-phase shift Butterworth), rectifying, and low-pass filtering (6 Hz, 4th order, zero-phase shift Butterworth) the raw EMG signals. We normalized all future EMG linear envelopes by the maximum value for each muscle measured during any of the maximum voluntary contraction activities. Participants then performed a standing static calibration trial for model scaling and hip circumduction trials for the identification of functional hip-joint-center locations 75 . To familiarize participants with the anatomy of their triceps surae muscles, we showed them images of the gastrocnemius and soleus muscles and palpated both muscles. To allow participants to feel the difference between gastrocnemius and soleus activation and to evaluate the differences in EMG signals, participants performed two sets of 10 doubleleg calf raises. The participants were standing (hip and knee at 0°) during the first set to primarily activate the gastrocnemius and were seated (hip and knee at 90°) with 4.5 kg weights on each knee during the second set to primarily activate the soleus. The average reduction in gastrocnemius-to-soleus activation ratio during the seated trial, compared to the standing trial, was 68% ± 14% for the 10 individuals who completed the study and 11% for the individual who was excluded from the study for the inability to reduce the activation ratio by at least 15% during the seated plantarflexion.
After acclimating to walking on the treadmill at 1.25 ms −1 for 5 min, participants performed five six-minute walking trials: a baseline trial, three feedback trials, and a retention trial (Fig. 2b). During the baseline walking trial, they were instructed to walk naturally. We computed their baseline gastrocnemius-to-soleus activation ratio (Eq. 3) and baseline stance-phase-averaged gastrocnemius EMG linear envelope. We averaged these values over the final minute of the baseline trial and normalized real-time measurements by the average baseline values in subsequent trials. During the first feedback trial, we used a real-time bar plot to instruct participants to reduce their activation ratio by at least 15% relative to baseline while maintaining normal walking kinematics. We instructed them to explore different strategies during the first four minutes of the trial, to converge on their most successful strategy during the fifth minute, and to maintain this strategy after feedback was removed for the sixth minute. The verbal instructions that were read to participants to explain the biofeedback are included in the Supplementary Information (Section: "Explanation of biofeedback given to participants during the experiment"). During the second and third feedback trials, participants were instructed to continue reducing their activation ratio below the target line, which was set to either a 15% reduction from baseline or their average reduction from baseline during the previous trial, whichever was greater. The 15% target was chosen as a feasible goal from pilot testing, but the adaptive target-setting method allowed for increasingly challenging goals for participants who were successful during the previous trial. During these final two trials, participants were given additional feedback of their average gastrocnemius activity, represented by the color of the bar. We instructed them to reduce their activation ratio by reducing gastrocnemius activation, or to achieve a small, blue bar (Fig. 2b). The total duration of feedback (15 min over 18 min of walking) was chosen based on the 15-20 min training duration used during the initial sessions of other gait retraining studies 41,76 . During the retention trial, we removed the feedback but instructed participants to walk with the coordination pattern that they learned during the feedback trials. A one-minute break was provided between each trial. We analyzed the final 30 steps of the sixth minute of all trials. For analysis, we used the same EMG filtering and normalization processes as used for the real-time experiment, and we averaged muscle activity over the duration of the stance phase.
EMG-informed simulation of knee contact force. We simulated these experimental data to estimate joint kinematics, kinetics, and KCF using OpenSim 61 and the previously described custom static optimization implementation. For simulations, five gait cycles were selected from each of the baseline and retention trials that had the smallest absolute differences in activation ratio compared to the average value from the final 30 steps in the trial.
Changes in plantarflexor EMG were incorporated into the static optimization simulation by constraining the simulated gastrocnemius-to-soleus activation ratio to match the ratio measured with EMG. First, a 40 ms electromechanical delay 72 was added to the medial gastrocnemius, lateral gastrocnemius, and soleus EMG linear envelopes. The simulated medial gastrocnemius-to-soleus activation ratio was constrained to match the EMG ratio within 2% (Eq. 7). A similar constraint was enforced for the lateral gastrocnemius-to-soleus activation ratio.
To validate our simulations, we qualitatively compared simulated activations to EMG linear envelopes with the 40 ms electromechanical delay (Fig. 5). For this comparison, we computed weighted averages of simulated activation and EMG linear envelopes for the vasti (vastus medialis and vastus lateralis) and biarticular hamstrings (semitendinosus and biceps femoris long head). The maximum isometric muscle force from the musculoskeletal model was used as the weight for each muscle.
We used the static optimization solutions and the Joint Reaction Analysis tool in OpenSim to compute KCF and decompose it into its muscle and intersegmental force components. We found the intersegmental reaction force by computing KCF from a simulation with all degrees of freedom actuated by ideal torque actuators. To elucidate the contribution of different muscle groups to KCF (Figs. 1, 6c), we prescribed the static (7) EMG med gastroc (t k − 40 ms) EMG med gastroc (t k − 40 ms) + EMG soleus (t k − 40 ms) − a med gastroc (t k ) a med gastroc (t k ) + a med gastroc (t k ) < 0.02 www.nature.com/scientificreports/ optimization-based muscle forces to a subset of muscles, generated the remaining net joint moments with ideal torque actuators, and subtracted the intersegmental force from the computed KCF. We defined the following functional groups of knee-crossing muscles: the quadriceps (vastus medialis, vastus intermedius, vastus lateralis, and rectus femoris), the hamstrings (biceps femoris long and short heads, semitendinosus, semimembranosus, gracilis, and sartorius), the gastrocnemius (medial and lateral gastrocnemii), and the tensor fascia latae.
Statistics. All statistical analyses were performed in MATLAB (R2017b) unless otherwise noted. After testing data for normality using the Shapiro Wilk test 77 , we compared normally distributed data using a two-sided, paired t test and compared non-normally distributed data using a two-sided Wilcoxon signed rank test. For comparisons of activation ratio and muscle activation changes across the feedback and retention trials (four comparisons), we report p-values after controlling for the false detection rate using R (v3.5.3, R Foundation for Statistical Computing, Vienna, Austria) 78,79 to avoid an increased Type I error rate. We used a two-one-sided t tests procedure to test for equivalence using the tool provided by Lakens et al. 80 with equivalence bounds of one baseline standard deviation to compare the average ankle plantarflexion moment between the baseline and retention trials. We defined peaks in the knee and hip joint angle and moment curves as the maximum value during the first 50% of the stance phase and the minimum value during the final 50% of the stance phase. Due to the exploratory nature of evaluating changes in peak joint angles and moments between baseline and retention, we did not correct for multiple comparisons. Values reported are mean ± standard deviation, and α = 0.05.

Code availability
The static optimization code and passive muscle force calibration code are available at https:// github. com/ stanf ordnm bl/ Matla bStat icOpt imiza tion and https:// github. com/ stanf ordnm bl/ Passi veMus cleFo rceCa libra tion, respectively. The real-time EMG biofeedback code is specific to our motion capture system and will be shared upon reasonable request to the authors. All code and data requests should be addressed to S.D.U. www.nature.com/scientificreports/